C      SUBROUTINES REQUIRED ARE READM, PRINTM, STAND, RCOEF AND SLE.
C
C      ================================================================
       DIMENSION X(100,20),XM(100,10),D(100,3)
       DIMENSION A(20,20),B(20),C(20)
       ND=100
       MD=20
       MM=20
C      INPUT MATRIX HAS N ROWS AND M COLUMNS (N= NO. OF OBS.)
C      AND M=NO. OF VARIABLES
       READ(3,*)N,M
       READ(3,*)((X(I,J),J=1,M),I=1,N)
       WRITE(4,*)' THE INPUT MATRIX IS: '
       WRITE(4,666)((X(I,J),J=1,M),I=1,N)
  666  FORMAT(7F8.1)
C
C      STANDARDIZE INPUT MATRIS AND THEN PRINT
C
       DO 201 I=1,N
       DO 201 J=1,M
       XM(I,J)=X(I,J)
 201   CONTINUE
       CALL STAND(XM,N,M,ND,MD)
       WRITE(4,*)' STANDARDIZED INPUT MATRIX: '
       DO 690 I=1,N
       WRITE(4,667)I,(XM(I,J),J=1,M)
 667   FORMAT(10X,I4,3X,7F8.1)
 690   CONTINUE
       CALL RCOEF(XM,N,M,ND,MD,A,MM)
       WRITE(4,*)'CORRELATION MATRIX. VARIABLE 1 IS Y'
       DO 22 I=1,M
       WRITE(4,668)I,(A(I,J),J=1,M)
 668   FORMAT(10X,I3,3X,7F10.3)
  22   CONTINUE
C
C      SET UP AND SOLVE SIMULTANEOUS EQUATIONS
C
       DO 100 I=2,M
       C(I-1)=A(I,1)
       DO100 J=2,M
       A(I-1,J-1)=A(I,J)
 100   CONTINUE
C
C      SOLVE SLE
C
       CALL SLE(A,C,M-1,MM,1.0E-08)
C
C      CALCULATE PARTIAL REGRESSION COEFFICEINTS
C
       DO 101I=1,M
       A(1,I)=0.0
       A(2,I)=0.0
       DO 101 J=1,N
       A(1,I)=A(1,I)+X(J,I)
       A(2,I)=A(2,I)+X(J,I)**2
 101   CONTINUE
       AA=N
       AB=N-1
       AC=SQRT((A(2,1)-A(1,1)*A(1,1)/AA)/AB)
       B(1)=A(1,1)/AA
       DO 102 I=2,M
       B(I)=C(I-1)*AC/SQRT((A(2,I)-A(1,I)*A(1,I)/AA)/AB)
       B(1)=B(1)-B(I)*A(1,I)/AA
 102   CONTINUE
C
C      CALCULATE ESTIMATED VALUE AND DEVIATION FOR EACH OBSERVATION
C
       DO 103 I=1,N
       D(I,1)=X(I,1)
       D(I,2)=B(1)
       DO 104 J=2,M
       D(I,2)=D(I,2)+B(J)*X(I,J)
 104   CONTINUE
       D(I,3)=D(I,1)-D(I,2)
 103   CONTINUE
       WRITE(4,692)
 692   FORMAT(10X,'COL 1 = Y,  COL 2 =ESTIMATED Y,    COL3 = DEVIATION')
       DO 693 K=1,N
       WRITE(4,669)D(K,1),D(K,2),D(K,3)
 669   FORMAT(12X,F6.2,10X,F10.5,10X,F10.5)
 693   CONTINUE
C
C      PRINT PARTIAL REGRESSION COEFFICIENTS
C
       WRITE(4,*)' REGRESSION COEFFICIENT   1 = CONST. TERM '
       WRITE(4,660)(B(K),K=1,M)
  660  FORMAT(7(E11.4,1X))
C
C      PRINT STANDARD PARTIAL REGRESSION COEFFICIENTS
C
       MMM=M-1
       WRITE(4,697)
  697  FORMAT(' STANDARD PARTIAL REGERSSION COEFFICIENTS ')
       WRITE(4,661)(C(K),K=1,MMM)
  661  FORMAT(6E10.4)
C
C      CALCULATE ERROR MEASURES
C
       SY=0.0
       SYY=0.0
       SYC=0.0
       SYYC=0.0
       DO 105 I=1,N
       SY=SY+D(I,1)
       SYY=SYY+D(I,1)**2
       SYC=SYC+D(I,2)
       SYYC=SYYC+D(I,2)**2
 105   CONTINUE
       SST=SYY-SY*SY/FLOAT(N)
       SSR=SYYC-SYC*SYC/FLOAT(N)
       SSD=SST-SSR
       NDF1=M-1
       AMSR =SSR/FLOAT(NDF1)
       NDF2=N-M
       AMSD=SSD/FLOAT(NDF2)
       R2=SSR/SST
       R=SQRT(R2)
       F=AMSR/AMSD
       NDF3=N-1
C
C      PRINT ERROR HEADINGS AND MEASURES
C
       WRITE(4,680)
 680   FORMAT('SOURCE OF',13X,'SUM OF   DEGREES OF   MEAN')
       WRITE (4,681)
 681   FORMAT('VARIATION', 13X,'SQUARES  FREEDOM    SQUARES     F-TEST')
       WRITE(4,683)SSR,NDF1,AMSR,F
 683   FORMAT(3X,'REGRESSION',5X,E11.4,I6,4X,E11.4,1X,E11.4)
       WRITE(4,684)SSD,NDF2,AMSD
 684   FORMAT(3X,'DEVIATION',6X,E11.4,I6,4X,E11.4)
       WRITE(4,685)SST,NDF3
 685   FORMAT(3X,'TOTAL VARIATION',1X,E11.4,I8,//)
       WRITE(4,686)R2,R
 686   FORMAT(3X,'GOODNESS OF FIT =',F10.4,5X,'CORR. COEFF. = ',F10.4)
       CALL EXIT
       END
       SUBROUTINE SLE(A,B,N,N1,ZERO)
       DIMENSION A(N1,N1),B(N1)
       DO 100 I=1,N
       DIV=A(I,I)
       IF (ABS(DIV)-ZERO)99,99,1
 1     DO 101 J=1,N
       A(I,J)=A(I,J)/DIV
 101   CONTINUE
       B(I)=B(I)/DIV
       DO 102 J=1,N
       IF(I-J)2,102,2
 2     RATIO = A(J,I)
       DO 103 K=I,N
       A(J,K)=A(J,K)-RATIO*A(I,K)
 103   CONTINUE
       B(J)=B(J)-RATIO*B(I)
 102   CONTINUE
 100   CONTINUE
       RETURN
 99    CALL EXIT
       END
C
C      SUBROUTINE TO CALCULATE THE MATRIX OF CORRELATIONS
C      BETWEEN COLUMNS OF DATA MATRIX X
C
       SUBROUTINE RCOEF(X,N,M,N1,M1,A,M2)
       DIMENSION X(N1,M1),A(M2,M2)
       AN=N
C
C      CALCULATE CORRELATION COEFICIENT BETWEEN COLULMNS I AND J
C
       DO 100 I=1,M
       DO 100 J=I,M
C
C      ZERO SUMS
C
       SX1=0.0
       SX2=0.0
       SX1X1=0.0
       SX2X2=0.0
       SX1X2=0.0
C
C      CALCUALTE SUMS, SUMS OF SQUARES AND SUM OF CROSS PROCUCT
C      OF COLUMN I AND J
C
       DO 101 K=1,N
       SX1=SX1+X(K,I)
       SX2=SX2+X(K,J)
       SX1X1=SX1X1+X(K,I)**2
       SX2X2=SX2X2+X(K,J)**2
       SX1X2=SX1X2+X(K,I)*X(K,J)
 101   CONTINUE
C
C      CALCULATE CORRELATION COEFFICIENT AND STORE IN MATRIX A
C
       RR1=(SX1X2-SX1*SX2/AN)
       RR2=SQRT((SX1X1-SX1*SX1/AN)*(SX2X2-SX2*SX2/AN))
       R=RR1/RR2
       A(I,J)=R
       A(J,I)=R
 100   CONTINUE
       RETURN
       END
C
C      SUBROUTINE TO STANDARDIZE THE COLUMNS OF A DATA MATRIX
C
       SUBROUTINE STAND(X,N,M,N1,M1)
       DIMENSION X(N1,M1)
       WRITE(*,*)'SUCCESSFUL ENTER'
C
C      STANDARDIZE EACH COLUMN OF THE MATRIX
C
       DO 100 I=1,M
C
C      CALCULATE MEAN AND STANDARE DEVIATION OF COLUMN
C
       SX=0.0
       SXX=0.0
       DO101 J=1,N
       SX=SX+X(J,I)
       SXX=SXX+X(J,I)**2
 101   CONTINUE
       XM=SX/FLOAT(N)
       SD=SQRT((SXX-SX*SX/FLOAT(N))/FLOAT(N-1))
C
C      SUBTRACT MEAN FROM EACH ELEMENT IN COLUMN, THEN
C      DIVIDE RESULT BY THE STANDARD DEVIATION.
C
       DO 102 J=1,N
       X(J,I)=(X(J,I)-XM)/SD
 102   CONTINUE
 100   CONTINUE
       RETURN
       END
